Brownian motion and Diffusion
Contents
Brownian motion and Diffusion#
import numpy as np
import scipy
from numpy.random import normal, choice, uniform
import ipywidgets as widgets
import matplotlib.pyplot as plt
plt.style.use('fast')
%matplotlib inline
Mean square displacement (MSD) of a random walker#
After time n number of steps (or time t) how far has random walker moved from the origin?
We quantify this by computing Mean Square Displacement (MSD). Note that the mean is computed over N number of simulated trajectories (ensemble average). Invoking central limit theorem, or simply realizing that off diagonal terms drop off we end up with expected 1/2 saling with steps.
Mean square displacement (MSD) and Diffusion#
Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement of particle in d dimensions to time \((t = n \Delta t)\).
Consider recording n number of positions for a particle in solution starting from origin \(r_0=0\)
Repeating experiment N number of times we can compute ensemble average
Where we decomposed displacement in d independent dimensions.
We next write number of steps in terms of time increments \(n=t/\delta t\)
Grouping constants together we defined the diffusion coefficient:
We end up with a general expression for a mean square displacement as a function of time. Any motion which adheres to this scaling will be called “diffusive”
MSD of a 1D random walk#
def rw_1d(T, N):
'''
L: trajectory length
N: Number of trajecotries
returns np.array with shape (T, N)
'''
# Create random walks
r = choice([-1,1], size=(T, N))
#Accumulate position
rw = r.cumsum(axis=0)
#Set initial position
rw[0,:]=0
return rw
T, N = 2000, 1000
rw = rw_1d(T, N)
t = np.arange(T)
R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time
msd = np.mean(R2, axis=1) # Notice we average over N
plt.loglog(t, np.sqrt(msd), lw=3)
plt.loglog(t, np.sqrt(t), '--')
plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);
MSD of a 2D random walk#
def rw_2d(T, N):
'''2d random walk function:
T: trajectory length
N: Number of trajecotry
returns np.array with shape (T, N)
'''
verteces = np.array([(1, 0),
(0, 1),
(-1, 0),
(0, -1)])
rw = verteces[choice([0,1,2,3], size=(T, N))]
rw[0, :, :] = 0
return rw.cumsum(axis=0)
traj = rw_2d(T=10000, N=1000)
#Simulate 2D random walk
T, N = 10000, 1000
traj = rw_2d(T, N)
#Compute RSD
dx = (traj[:, :, 0]- traj[0, :, 0])
dy = (traj[:, :, 1]- traj[0, :, 1])
R2 = np.mean(dx**2 + dy**2, axis = 1) # notice how we averaging over N
fig, ax = plt.subplots(nrows=2, figsize=(10,10))
t = np.arange(T) # time axis
ax[1].loglog(t, np.sqrt(R2), lw=3, alpha=0.5);
ax[1].loglog(t, np.sqrt(t), '--');
ax[0].set_title('2D random walker',fontsize=15)
ax[0].plot(traj[:3000, :5, 0], traj[:3000, :5, 1]);
ax[0].set_xlabel('X')
ax[1].set_xlabel('Y')
ax[1].set_xlabel('Number of steps, n',fontsize=15)
ax[1].set_ylabel(r'$MSD(n)$',fontsize=15);
Brownian motion#
The Brownian motion describes the movement of a particle suspended in a fluid resulting from random collisions with the quick molecules in the fluid (diffusion).
A small particle undergoes large number of molecular collisions when going from one step to another or after \(dt\) time. As a result each displacement over time \(dt\) can be viewed as a sum of ranomd collisions which can be approximated by a normal distribution via Central Limit Theorem.
Thus to simulate brownian motion we draw ranom displacements from normally distribution.
We assume we have started at position \(\mu=0\) and our variance is given by \(\sigma^2=2Dt\) Where D is the diffusion coefficnets which is related to parameters of discree random walk as shown in the lecture.
In the last step we re-wrote brownian motion equation in a convenient way by shifting normally distributed radnom variable by \(\mu\) and scaling it by \(\sigma\)
def brown(T, N, dt=1, D=1):
"""
Creates 3D brownian path given:
time T
N=1 trajecotires
dt=1 timestep
D=1 diffusion coeff
returns np.array with shape (N, T, 3)
"""
n = int(T/dt) # how many points to sample
dR = np.sqrt(2*D*dt) * np.random.randn(N, n, 3) # 3D position of brownian particle
R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
return R
Below we proceed to simulate continuous random walk in 1D-3D
We will plot trajecotires and distributions of brownian particle using interactive plotting via ipywidgets and holoviews/plotly interface.
R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
@widgets.interact(t=(10,3000-1))
def brownian_plot(t=10):
fig, ax = plt.subplots(ncols=2)
ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
ax[1].hist(R[:, 10, 0], density=True, color='red')
ax[1].hist(R[:, t, 0], density=True)
ax[1].set_ylim([0,0.1])
ax[0].set_ylim([-200, 200])
ax[0].set_xlim([-200, 200])
fig.tight_layout()
import holoviews as hv
hv.extension('plotly')
plots = []
for i in range(10):
plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
plots.append(plot)
hv.Overlay(plots)
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]
xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])
hv.Overlay(rw_curve) << ydist << xdist
Diffusion Equation#
The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\)
Diffusion equation:
Formulated empirically as Fick’s laws
This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour
This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.
Can solve numerically by writing derivatives as finite differences!
Can also simulate via random walk!
Diffusion coefficient \(D\), Units \([L^2]/[T]\)
Important special case solution (here written in 1d):
where \(\sigma(t)=\sqrt{2{D}t}\)
density remains Gaussian
Gaussian becomes wider with time
check that this is indeed a solution by plugging into the diffusion equation!
def sigma(t, D = 1):
return np.sqrt(2*D*t)
def gaussian(x, t):
return 1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
@widgets.interact(t=(1,100,1))
def diffusion(t=0.001):
R = brown(T=101, N=1000)
x = np.linspace(-20, 20, 100)
plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
plt.legend()
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
plt.xlim([-25, 25])
References#
The mighty little books
More in depth
On the applied side